
#include <gsl_ext_oct/OctaveArgParser.h>

#include <octave/oct.h>

const string ERROR_PREFIX = "GSLpp_testfunction_data: ";

bool checkArgs(const octave_value_list &args)
{
  const int nargin = args.length();
  
  if(nargin != 5)
  {
    error((ERROR_PREFIX + "5 parameters required.").c_str());
    return false;
  }
  
  if(!args(0).is_string())
  {
    error((ERROR_PREFIX + "The first parameter (function name) must be a string.").c_str());
    return false;
  }
  
  if(!args(1).is_map())
  {
    error((ERROR_PREFIX + "The second parameter (function parameters) must be a struct.").c_str());
    return false;
  }
  
  if(!args(2).is_real_type())
  {
    error((ERROR_PREFIX + "The third parameter (number of mesh points) must be an integer.").c_str());
    return false;
  }
  
  if(!args(3).is_matrix_type())
  {
    error((ERROR_PREFIX + "The fourth parameter (x range) must be a vector.").c_str());
    return false;
  }
  
  if(!args(4).is_matrix_type())
  {
    error((ERROR_PREFIX + "The fifth parameter (y range) must be a vector.").c_str());
    return false;
  }
  
  return true;
}

DEFUN_DLD(GSLpp_testfunction_data, args, nargout, " -*- texinfo -*-\n\
@deftypefn {Loadable Function} {[@var{x},@var{y},@var{fxy}] =} GSLPpp_testfunction_data (@var{f},@var{fp},@var{n},@var{xrange},@var{yrange})\n\
Evaluates the given function @var{f} in an \n\
@var{n}x@var{n} grid. Additional function parameters are \n\
specified with the struct @var{fp}. The ranges of x- and y-\n\
coordinates are specified with two-value vectors \n\
@var{xrange} and @var{yrange}, respectively.\n\
@end deftypefn")
{
  octave_value_list retval;
  
  if(!checkArgs(args))
    return retval;
  
  string funcName = args(0).string_value();
  Octave_map funcParams = args(1).map_value();
  int nPoints = args(2).scalar_value();
  ColumnVector xRange = ColumnVector(args(3).vector_value());
  ColumnVector yRange = ColumnVector(args(4).vector_value());
  
  Octave_map funcEvalParams;
  Matrix X, Y, FXY;
  
  funcParams.assign("n", 2);
  
  funcEvalParams.assign("f", "sym");
  funcEvalParams.assign("g", "fd");
  funcEvalParams.assign("H", "fd");
  
  GSLFunctionWrapper *f;
  gsl_multimin_function *f_gsl;
  gsl_ext_testfunction_params tfGSLParams;
	
  try {
    f = OctaveArgParser::getFunctionInstance(funcName,
                                             funcParams,
                                             funcEvalParams,
                                             GSLFunctionWrapper::F,
                                             &tfGSLParams);
    f_gsl = dynamic_cast< GSLFunctionWrapper_f * >(f)->
            operator gsl_multimin_function*();
  }
  catch(exception &e) {
    error((ERROR_PREFIX + "Unable to initialize the given function.").c_str());
    error(("Reason: " + string(e.what())).c_str());
    return retval;
  }
  
  X = Matrix(nPoints, nPoints);
  Y = Matrix(nPoints, nPoints);
  FXY = Matrix(nPoints, nPoints);
  gsl_vector *x = gsl_vector_alloc(2);
  
  for(int i = 0; i < nPoints; i++)
  {
    for(int j = 0; j < nPoints; j++)
    {
      X(i, j) = xRange(0) + j * (xRange(1) - xRange(0)) / (nPoints - 1);
      Y(i, j) = yRange(0) + i * (yRange(1) - yRange(0)) / (nPoints - 1);
      gsl_vector_set(x, 0, X(i, j));
      gsl_vector_set(x, 1, Y(i, j));
      FXY(i, j) = GSL_MULTIMIN_FN_EVAL(f_gsl, x);
    }
  }
  
  gsl_vector_free(x);
  delete f;
  
  retval(0) = X;
  retval(1) = Y;
  retval(2) = FXY;
  
  return retval;
}
